#include "VGQuantileBounds.h"
#include <math.h>
#include <boost/math/special_functions/gamma.hpp>
//#include <boost/math/constants/constants.hpp>

VGQuantileBounds::VGQuantileBounds(double mu, double lambda, double alpha, double beta) 
	: _mu(mu), _lambda(lambda), _alpha(alpha), _beta(beta), _asympExpansion(mu, lambda, alpha, beta, 0)
{
	using namespace boost::math;

	_lowerCutoff = pow(2,-lambda)*pow(alpha,-lambda)*exp((alpha - beta)*(-1 + lambda)*pow(alpha + beta,-1))*pow(tgamma(lambda),-1)*
		pow(pow(alpha,2) - pow(beta,2),lambda);

	_upperCutOff = pow(2,-lambda)*pow(alpha,-lambda)*exp(-1.0)*pow(tgamma(lambda),-1)*
		(exp(1.0)*tgamma(lambda)*pow(2,lambda)*pow(alpha,lambda) - 
		exp(lambda)*pow(pow(alpha,2) - pow(beta,2),lambda));
}


VGQuantileBounds::~VGQuantileBounds(void)
{
}

std::pair<double, double> VGQuantileBounds::operator()( double u ) const
{
	return std::make_pair (LowerBound(u), UpperBound(u));
}

double VGQuantileBounds::UpperBound( double u ) const
{
	if(u >= _upperCutOff)
		return _asympExpansion.RightTail(u) + 1;
	
	return _asympExpansion.RightTail(_upperCutOff);
}

double VGQuantileBounds::LowerBound( double u ) const
{
	if(u <= _lowerCutoff)
		return _asympExpansion.LeftTail(u) - 1;

	return _asympExpansion.LeftTail(_lowerCutoff);
}
